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We present an exact Monte Carlo algorithm designed to sample theories where the energy is 
a sum of many couplings of decreasing strength. Our algorithm, simplified from that of L. Lin et 
al. JO, avoids the computation of almost all non-leading terms. We illustrate its use by simulating 
5(7(2) lattice gauge theory with a 5- loop action, and discuss further applications to full QCD. 



I. INTRODUCTION 
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When sampling by Monte Carlo the partition function Z = JJ\dU e H ^ u ^\ the most common algorithm is 
that of Metropolis g: at each step, starting from the current configuration {U}, a candidate configuration {U 1 } is 
bJ[). proposed, and it is accepted with probability 
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min(l,e-( H « u '»- H « u »)) . (1) 



This acceptance test is realized by comparing the right-hand side of (|lj) to a random number uniformly distributed in 
[0, 1]. This seems like a waste of information: why compute H({U'}) exactly, then compare it with a random number? 
It should be sufficient to estimate it. Indeed, this logical proposition has been studied several times Two 
\Q , difficulties have been identified, both caused by the non-linear relationship between the energy H and the probability 
oc e~ H : (i) what is needed is an unbiased estimate of e~ H , which must be obtained from unbiased estimate(s) of H; 
(ii) to be interpreted as a probability, the noisy estimator of P aC c must be bounded, and in particular stay positive. 
Difficulty (i) was overcome in Q, which however showed that violations of (ii) caused intolerable systematic errors 
unless the amount of noise in the estimate of H was minuscule. Difficulty (ii) was overcome in Q, which showed that 
exact results could be obtained even in the presence of a large amount of noise in the estimate of H. Ref. Q, however, 
introduces an infinite number of auxiliary variables, so that it may take an infinite amount of work to bring these 
auxiliary variables to equilibrium. And tests of the method are performed on a toy model with 5 degrees of freedom 
'"" p 1 1 only, whose relevance may be questioned. Here, we simplify the method of [Q , by introducing only 1 auxiliary variable 
O I per term in H. Moreover, we separate H into a leading part to be calculated exactly, and a sum of small correction 
terms, which we treat stochastically. This separation is essential: because stochastic estimates are used for correction 
terms only, large amounts of noise can be tolerated. As a consequence, our algorithm is a very efficient approach to 
the simulation of complicated Hamiltonians. 
■ Consider a generic Hamiltonian of the type 

b ■ 

H = Y J C k W k (2) 

where as k increases, \c k \ decreases and the successive terms W k typically become less and less local. For instance, 
in a spin model {o^}, Wq would be the nearest-neighbour interaction X)<ij> Wi would represent next-nearest- 

neighbour interactions, etc... Here, we will illustrate our method for lattice gauge theory. In that context, W k are 
the traces of Wilson loops of increasing size: Wq = J2 X u u ^ r EL ^ around elementary plaquettes, Wj.,2,3 correspond 
to different geometries of 6-link loops, etc... It is often the case that one would like to study a Hamiltonian of 
type (g) resulting from an expansion, be it perturbative ]^], non-perturbative ||, or based on the fixed point of a 
renormalization group transformation [Q. In all these situations, the expansion is truncated to a maximal order m 
dictated by technical reasons. As k increases in (^|), the number of geometrically equivalent terms grouped into W k 
increases exponentially: in a spin model on a hypercubic lattice in d dimensions, each spin has 2d nearest-neighbours 
(these interactions are grouped into Wo); d(d— l)/2 next-nearest neighbours (grouped into W\), d(d— l)(d — 2)/6 3rd- 
neighbours, etc... This combinatoric explosion normally makes the simulation of extended Hamiltonians prohibitively 
expensive. This is the reason for a truncation to very low m, often taken to be 1 or 2. However, in most cases the 
couplings c k in (j^) decrease exponentially with k, so that the overall Hamiltonian is dominated by Wq, with small 
corrections. In lattice field theory, this is actually required if the Hamiltonian is to make sense and tend to a local 
operator as the continuum limit is approached. By making use of stochastic methods to estimate the correction terms 
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Wk, k > 1, we aim at postponing the combinatoric explosion of the simulation costs incurred when including higher 
terms Wk- This opens the possibility of studying numerically much more complicated Hamiltonians including higher- 
order correction terms. In lattice field theory, these correction terms are crucial to suppress discretization errors, and 
form the building blocks of so-called "improvement" strategies. Also the inclusion of higher-order terms can be very 
useful in the approaches to the fermion determinant simulations involving the loop expansion 

We present our method in Section II, and illustrate it in Section III with simulations of a 5-loop perturbatively 
improved action for SU (2) lattice theory. We conclude with prospective applications of our method, in particular for 
dynamical fermion QCD simulations. 



II. NOISY MONTE CARLO: THE METHOD 

Given the Hamiltonian (^) let us suppose that the terms CkWk are nonpositive starting from k = 1: 

k > 1 : c k W k (U) < W (3) 

This can be easily arranged by adding to each term of the Hamiltonian a nonessential constant. Here U are the 
fields of the model under consideration. The key idea of the method is to estimate the contribution of the terms 
Wk(U), k > 1 stochastically by introducing auxiliary fields. This will lead to a significant reduction of computational 
effort if the coefficients Cfe, k > 1 are small enough. In all cases, the algorithm remains exact. 

We introduce auxiliary fields a k , k > 1 (associated with the terms Wk), which can take two values: and 1. Using 
the identity 

a + b= ^2 [a*5 a ,Q + b*da,i] (4) 

cr=0,l 

we represent the probability e~ H in the form 

e- H = P [U]*P 1 [U,a] (5) 

where 

m 

P [U]^e-^o(U) . P 1 [U,a} = m2l^,o+S ak A^ ChWkiU) - 1 )] (6) 

k = l <Tfc=0,l 

The r.h.s. of (||) can be interpreted as the joint probability distribution for the original fields of the model and the 
new a fields. Due to the inequalities (||) this distribution is well defined: Pi[U, a] > V{Z7, a}, and the probabilities 
for Gk to take value or 1 when the U fields are fixed lie in the interval [0,1]: 

This means that our algorithm has no probability bound violations, which plagued previous attempts to construct an 
efficient stochastic algorithm J|,|]]. 

One can easily see why the introduction of auxiliary a fields can be useful. Starting from the current {U\,a} 
configuration, a candidate configuration {U2} distributed with the weight Po[C^2] is proposed, and accepted with 
probability 

p • A. Pi[U 2 ,^h . /, -,-r e-^W-l x 

Pace = mm 1 . — r = mm 1 , ,,. , ( 8 ) 

V 'Pi Ui,o- / V ' J-J- e - c k w k(Ui) _ \) K ' 



k : (Tk = l 



Since the terms CkWk(U) contribute in P acc only if 07. = 1, the amount of computational work is greatly reduced if 
the configurations with Ok — are dominating. That is certainly the case when the absolute values of the coupling 
coefficients \ck\ are small: the probabilities for <ik to be unity, averaged over {U} configurations, are negligible then. 
Indeed, to leading order in Ck the average probability Pa k =i from eq. (j7|) can be written as 

(Pa*=i) « -c k {W k {U)) wO if c k ^0 (9) 

Expression ([)]) also suggests that one should try to make |(Wfe)| as small as possible, using the freedom one has to 
shift Wk by a constant. This goal should remain compatible however with inequalities (^); otherwise, probability 
bound violations will appear for p ak and P acc in eqs.(0j|). 
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Actually, the violation of conditions (|J) is not completely forbidden. As it was pointed in Ref. Q, one can address the 
problem of the lower probability-bound violations by redefining the measure. If the distribution Pi[U, a] in (^) can be 

negative for some configurations {U, er}, one can effectively simulate with the probability distribution Pq[U] * Pi [U, <r] 

instead and include the sign sgn{P\) into the observable expectation value: 

{ } (sgn(Pi))\\ [ ' 

where by ()|| we denote the averages with respect to distribution Po[f] * Pi[U, cr] ■ Sometimes the admission of very 
rare sign violations can substantially decrease the probability p ak =i- However one should be very careful in using this 
trick: as the volume of the system increases, one needs an exponential growth of statistics to estimate (sgn(Pi))n 
within the same accuracy. In the following we shall always assume fulfillment of the inequalities ([5|). 

After updating the U fields one should also update the a fields to preserve ergodicity. This requires the calculation 
of probabilities (0). At this point the reader might say: "OK, one saves computational effort by not calculating some 
terms in expression (^) while estimating P acc . Nevertheless, one must calculate these terms when updating the a 
fields! So does one gain anything in the end?" The answer is 'y es ' f° r the following two reasons: 

• The terms Wu for which it is reasonable to use the stochastic estimation usually couple many degrees of freedom 
(this is due to the usual nonlocality of weakly coupled terms, which serve as corrections to more local leading 
terms in the Hamiltonian) . If one uses usual local algorithms (without introducing stochastic a variables), one 
should estimate the term Wk each time one updates a degree of freedom which it couples. Contrary to that, if 
one uses Noisy Monte Carlo, the probabilities (R) should be calculated only once per a update. 

• The variables can be refreshed infrequently, the more so as the associated coupling c& gets smaller. It will 
be demonstrated in the next Section on a particular example. This slow dynamics of the auxiliary a fields does 
not imply slow dynamics of the physically relevant U fields. 

Up to now we were quite generic, showing that the Noisy Monte Carlo (NMC) method can be potentially very 
effective for the variety of theories where the energy (|2j) is a sum of couplings of decreasing strength. In the next 
section we illustrate these ideas on a particular example: a 5-loop perturbatively improved SU(2) Yang-Mills model. 



III. 5-LOOP SU(2) GAUGE THEORY 

We consider a 5-loop SU(2) gauge action in Ad: 

5 1 

— ' m^r 



<=i m i n i 



(11) 



where the indices (m^rij) = (1, 1), (2, 2), (1, 2), (1, 3), (3, 3) for i = 1,...,5 denote the planar, fundamental loops of 
size m x n 



S ,i 



]T ( -2 * sgn ( Ci ) - ^ ( 



(12) 



X fJL 



The Gibbs factor is exp(— ^S). Note that in eq.([12|) we have arranged the constant term —2 * sgn (c,) to ensure the 
condition (^|) for elementary action terms corresponding to each loop: 
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-sgn ( Ci ) - — 



< 



(13) 



Using the results of H one can construct a one-parameter set of actions which have no 0(a 2 ) and C(a 4 ) corrections 



ci =(19 - 55 cb)/9, c 2 = (1 - 64 c 5 )/9 
c 3 =(-64 + 640 c 5 )/45, c 4 = l/5-2c 5 



(14) 



Here we take C5 = 1/20 (the same action was used in the context of improved cooling JTo]|). 
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Following the ideas of section II, we estimate the contribution of all loops except the plaquette stochastically. For 
each loop I = {fj,, v, x} of sort 2 < i < 5 we introduce the auxiliary variable oi (/) = 0, 1; and rewrite the contribution 
of this loop to Gibbs factor in the form 

e -f *.„,„.. = £ [5ff<(0i0 + j CTj(0)l(e -#^„_ 1)] (15) 

<Ti(0 = 0,l 

The resulting distribution of {U, a} fields is used for the generation of independent {U} configurations. We shall say 
that for a given {a} configuration the loop {I, i} is "active" if cr^ (Z) = 1. 

Let us describe the updating procedure in the {[/, a} configuration space. Consider first the local updating of the 
gauge fields U when the a fields are fixed. The proposal value U£ e ™ at a given link {x, fi} is generated by heatbath 
with respect to the measure 

P [U] cx exp(-£cxSi,i[U]) (16) 



Pace = minfl, - 1 . (17) 



where Si,i is the plaquette action (see (|12|)) ; and then accepted with probability 

1C\\ 



where 



Pi[U,a}= [] (e~i S -' [U] -l) (18) 

(/,«) 9 {x,fi} 

(Ti(l) = 1 



Only active loops which contain the given link contribute to the expression in the r.h.s. of (|lj 

After each Ni updates of fields U on the entire lattice we update the a fields of sort i. For each loop I we assign 
the values 0,1 to the variable with the following probabilities 

Pa z (i)=o = exp(-5i,i[C/]) ; Pa % (i)=i = 1 - exp(-Si t i[U]) (19) 



Due to the absence of interaction between different a variables, the probabilities (19) depend only on the gauge 
configuration, so that a variables can be updated independently. 

In our simulations we have measured the average values of <Ji, 2 < i < 5 which are listed in Table I. They are quite 
small, and very close to the perturbative estimate (^). This shows that one can avoid the computation of almost all 
of the extended "staples" in the U update. 



loop 


1x2 


1x3 


2x2 


3x3 




0.0753 


0.0199 


0.0202 


0.0018 



Table I: Average value of en field for each loop of sort i. 

Performing numerical simulations for the 5-loop model (|TT| , [l4l ) with auxiliary a fields, we were mainly interested in 
the efficiency of our new NMC algorithm. One can estimate the efficiency of NMC by comparing it with the updating 
procedures which are commonly used now for the simulation of multiloop actions like (|l l|) . In the following we label 
these usually applied techniques with the collective name "Usual Monte Carlo" (UMC), to contrast with NMC. 

We compare the computer times needed to get the same results with NMC and UMC algorithms as follows. First, 
we make an analytic estimation of the total computational cost of one update of the U fields for both algorithms in 
units of matrix (link) multiplications. Second, we extract from numerical simulations the integrated autocorrelation 
times for different observables, in units of U update. The computer time needed to estimate any given observable is 
proportional to the product of the computational cost per update and the autocorrelation time. 

For NMC the average computational cost of one update of the U fields on the entire lattice is equal to 

5 

t NMC = t Pj + 4Vj2 n *ta P le(i) * 7l mult (l) * fa) (20) 
i=2 
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where tfj is the cost for generating the proposal configuration with measure ( fHf ) (i.e. the update cost for the elementary 
plaquette action), AV is the number of links on the lattice, n stap ie(i) is the number of "staples" which the loops of 
sort i form for each link, n mu it{i) is the number of matrix multiplications needed to estimate the contribution of one 
staple of sort i, and the factor (a) accounts for the fact that one needs to calculate the contribution of active loops 
only. One can easily check that 

3 

n st aple{i) — ^ * P i * S * > n mult{i) = Pi (21) 

where Pi = 2(rm + rij) is the perimeter of loop i, and Sj is a symmetry factor: s,; = 1 for square loops and Sj = 2 for 
rectangular loops. Then wc have 

t NMC = t pl +&v j-p2 Si{(Ji) (22) 

»=2 

On the other hand, the computational cost of one update of the u\ fields on the entire lattice is given by 

t a% = (23) 

Here 6Us; is the number of loops of a given sort on the lattice, and the perimeter Pj of the loop appears again as the 
number of matrix multiplications n mu it{i) needed to calculate the probabilities ([l9|). Since we update the ct 2 ; fields 
only once per each Ni updates of the U fields, the total computational cost per U update for the NMC method is 



t NMC = t NMC + £ k = # + wj^PM^ + PM)) (24) 



i=2 



Let us note that the computational cost t^Jf for UMC of one U update is approximately equal to the r.h.s. of 
expression ( pi] ) in the limit A, ; — > oo and (crj) — > 1: 

t^ c = t^ + 6Vj2PS (25) 



Indeed, in the limit when all a are set equal to 1 and not updated we recover the usual algorithm (certainly one 
should correct the expressions ( |l7| , ^8|) for P acc in this case). 

Now we can compare the performance of our NMC algorithm with that of UMC. The naive gain in efficiency from 
using NMC does not depend on the observable measured, and is equal to the ratio between the computational costs 
© and ((§): 

f UMC fP'_i_fiT/V 5 P 2 S- 
naive _ L tot _ L U ' u v _ 2^ii=2 fj (n fi \ 

'gain — .NMC ,pl . C Tr^5 n / 1 , r> I \\ ^ ' 

hot t P v + 6V 5J i=2 PiSiij^ + Pi{(Ti}) 

Now one should also take into account the increase of autocorrelation times coming from the introduction of auxiliary 
variables a in the NMC algorithm, so that the real gain is 

UMC 

real _ naive int (07) 
'gain — 'gain * NMC \ > 

'int 

where t^ ic and t^ ic are integrated autocorrelation times for UMC and NMC respectively. Note that T^f 1 is a 
function of the updating frequencies 1/JVf of the a fields. Like Tint, the ratio ( p7| ) will also depend on the observable. 

In Table II we present the autocorrelation times for averaged traces of 6 different loops in units of U updates. In 
the first column we show the results for the Usual Monte Carlo, and in other columns for the NMC algorithm with 
different frequencies of a updates for 1x2, 1x3, 2x2, 3x3 loops. In the last row we present the naive gain (|2^). 

Let us make one useful remark. It is not necessary to keep the same updating frequencies 1/ATj for all sorts i of loops. 
Actually it is even impractical. The computational cost of U update coming from the loop of sort i is proportional to 
the average value of Cj, which is in turn proportional to the coupling (|^). As the coupling decreases, we should expect 
a reduction of the computational effort for the corresponding terms in the action. That is not the case for the cost 
of <7j update: it does not depend on the coupling and even increases with the nonlocality of the action term (factor 
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Pi in expression (|23|)). In order for the work in the a and in the U updates coming from loops of sort i to remain 
comparable, one should keep the updating frequencies l/Ni proportional to (oj): 

]Jr ~ Pi(*i) (28) 

Due to the small influence of weakly coupled terms on the dynamics of the system, one can expect only insignificant 
changes in the autocorrelation behavior as N{ increases. These considerations are distinctly demonstrated in Table 
II, where in two columns we present the results for updating frequencies of a fields varying in accordance with . 



number of 


UMC 


1 


5 


5 for 1x2 


10 


10 for 1x2 


20 


30 


40 


50 


U updates per 


no 


for 


for 


15 for 1x3,2x2 


for 


30 for 1x3,2x2 


for 


for 


for 


for 


1 a update 


a 


all 


all 


105 for 3x3 


all 


210 for 3x3 


all 


all 


all 


all 


T^t (1x1) 


0.7(1) 


1.9(2) 


2.3(1) 


2.5(2) 


3.1(2) 


3.2(2) 


4.3(4) 


4.5(4) 


3.8(2) 


5.1(4) 


Tint (1x2) 


0.8(1) 


2.6(3) 


2.8(2) 


3.2(2) 


4.3(4) 


3.9(3) 


5.2(4) 


5.6(5) 


5.7(4) 


7.7(8) 


T int (1x3) 


0.8(1) 


2.7(3) 


2.8(2) 


3.2(2) 


4.3(4) 


3.9(3) 


5.1(4) 


5.4(5) 


5.4(4) 


7.4(8) 


Tint (2x2) 


1.0(1) 


3.4(5) 


3.3(3) 


3.7(3) 


4.7(4) 


4.7(3) 


5.3(4) 


5.8(6) 


5.7(4) 


7.5(8) 


Tint (2x3) 


1.4(3) 


4.2(7) 


3.8(3) 


4.2(3) 


5.5(5) 


5.7(4) 


5.7(5) 


6.2(6) 


6.3(5) 


7.7(8) 


Tint (3x3) 


1.8(4) 


5.0(8) 


4.5(5) 


5.0(5) 


5.8(6) 


6.4(5) 


5.9(5) 


6.3(6) 


6.3(5) 


7.5(8) 


r naive 
qain 


1 


7.2 


14.8 


18.3 


17.9 


20.6 


18.3 


20.8 


21.1 


21.4 



Table II: Integrated autocorrelation times for average loop traces in units of U updates for UMC algorithm (first column) 
and for NMC algorithm with different frequencies of a updates for 1x2, 1x3, 2x2 and 3x3 loops (other columns). The last row 
presents the naive gain for the NMC algorithm (|2^). 

Table II gives an impressive demonstration of the benefits which come from using the NMC algorithm. The naive 
gain increases substantially as we decrease the frequencies of a updates, while the autocorrelation times grow rather 
slowly. That is particularly visible for the runs where the updating frequencies for a fields are adjusted as per eq. (p8|) . 
For such runs we can infer that the 'real gain' 0(4 — 6) in computer time ( f27|) for the observables measured is large 
enough for a convincing demonstration of the possible advantages coming from using the NMC algorithm. 

Let us make a conclusion for this section. We have applied our NMC algorithm for the 5doop model (|ll],|l4|). We 
have shown that with this algorithm a significant gain in efficiency is obtained in comparison with usual updating 
techniques. Finally we note that the action ( pd| ) is a relatively simple one, and one can expect a much greater gain 
for more complicated highly-improved actions with many nonlocal weakly coupled terms. 



IV. DISCUSSION 



Let us summarize our algorithm: 

a) Separate the Hamiltonian (or action) into a dominant term cqWq, to be calculated exactly, and correction terms 
Y^k=i c kWk, to be estimated stochastically. 

b) Shift the correction terms to guarantee CkWk < 0. 

c) Introduce auxiliary local variables o~k(l), through identity (^) (here I runs through all the elementary 'bonds' which 
form Wk, e.g. loops in gauge theory). 

d) Update the auxiliary variables Ok by heatbath. 

e) To update the original variables U, propose a new value U' sampled from the distribution oc e~ c ° w ° , and accept it 
with the Metropolis probability min(l, n k >i ;0 . k= i(e" CkWk(u,) - l)/(e" CkWk ( u ) - 1)). 

The essential advantage of our algorithm appears in step (e): only the terms Wk whose associated Uk is equal to 
1 need to be computed. Since on average, (o~k) goes to zero with Ck, the computation of almost all correction terms 
can be avoided. 

To avoid simply shifting the cost of the algorithm to step (d), we propose to refresh the variables o~k infrequently, 
the more so as the associated coupling \ck\ gets smaller. We have pointed out that this introduction of slow dynamics 
for the Uk does not enforce slow dynamics for the system, since Wk(l) will fluctuate regardless of the value of cr fc (Z). 
Our numerical study of Section III confirms this statement. 

Let us now speculate on possibilities to use our algorithm to simulate a Hamiltonian with a very large number of 
terms. A specific example we have in mind is the case of full QCD, where the measure is, for 2 flavors of Wilson 
quarks: 
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^e" s ^ u) det 2 {l- kM(U)) (29) 

where S g is the local gauge action, M(U) is a hopping matrix connecting nearest neighbours on a Ad hypercubic 
grid, and Z normalizes the distribution. The determinant can be turned into exp(Tr(Log(l — kM (U)))) , then the 
logarithm expanded around 1, giving the loop expansion of the measure above: 

l e - Sg (u)-2Y,ZA TrM W . (30) 

TrM(U) 1 can be represented as a sum over all closed non-backtracking loops of length I on the Ad hypercubic lattice. 
The number of types of contributing loops ni is bounded by (2d — 1)' — 7 l , because of the branching factor at each 
hop. Although this upper bound is not saturated, it is clear that the multiplicity of terms of a given length I grows 
exponentially: 

ni~Fi{l)a l ; a<7 (31) 

where F\(l) is a rational function of I, and a 1 is the leading exponential ascend of the number of loops of length I in 
the limit of large I. For this reason, it seems that sampling numerically the distribution (^0|) is a disastrous idea: the 
action contains an infinite number of terms, of exponentially growing multiplicity. Instead, other strategies are being 
used, based on the transformation of the determinant (^9|) into a Gaussian integral. 

Nevertheless, the coupling 4- decreases exponentially as I increases. Therefore, the auxiliary variables ai associated 
in our approach with various loops of length I will take value almost always. From (^|) one gets 

(m) F 2 {l)k l (3 l (32) 

where the exponentially growing factor (3 comes from the average trace of Dirac matrices along the loops of length Z, 
i*2(0 is again the rational function. 

If one arranges the updating frequencies for each loop i as per eq. (|28|) , one can expect that the average computer 
time needed for estimation of contribution of all loops of length I behaves at large I as 

ti ~ nil 2 (m) ~ F(l) * (af3n) 1 (33) 

(It was pointed above that one should not expect a significant growth of autocorrelation coming from slow dynamics 
for ai ) For k < K ca = -^3 the computational cost ti decreases exponentially with Z and the total computational cost 
of the algorithm t = ti converges to some finite value. Our rough estimation from fitting ni and average trace of 
Dirac matrices in the interval 4 < I < 12 gives a « 5.4; (3 « 1.4, and therefore K ca ~ 0.13. 

In the regime k < n ca we are in an interesting situation where the influence of very large loops is negligible because 
their associated coupling in the effective action is extremely small. Therefore, truncating the loop expansion above a 
certain order will introduce a statistically unobservable bias. Equivalently, one can freeze the associated a variables at 
the value zero, or update them with arbitrarily low frequency. In spite of this extremely (or infinitely) slow dynamic 
mode of the er's, the dynamics of the gauge fields are not affected. Note that the cost of our algorithm grows linearly 
with the volume V of the system. This is better than alternative approaches to the simulation of full QCD: Hybrid 
Monte Carlo (cost oc V 5 ^ 4 ) and MultiBoson (cost oc V(LogV) 2 ) p| . In addition the stepsize, or typical change at 
each update of a gauge link J7, does not seem restricted a priori for small quark mass, unlike in the two alternative 
approaches above. Unfortunately the possible high efficiency of our algorithm is counterweighted by its extreme 
programming complexity. 

A less speculative use of our algorithm for full QCD consists of truncatingthe loop expansion eq.(|30|) to some 
order l ma x, and to represent the higher orders with the MultiBoson approach |12|| . This strategy, called "UV-filtered 
MultiBoson" , has already been used successfully ||. However, in Ref. || the loop expansion is truncated to its lowest 
term / = 4, because the exact evaluation of larger loops is too time-consuming. With our stochastic approach, these 
larger loops can be estimated at low cost. We expect this composite strategy to be particularly efficient. 

Acknowledgments: We acknowledge communications with T. DeGrand, M. Ilgenfritz and U. Wenger. T.B. was 
supported by INTAS under grant 96-0370 and Russian Basic Research Fund under grant 99-01-00190. 
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